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Abstract 



The formation and evolution of nonlinear and turbulent dynamical structures in two-dimensional 
complex plasmas and fluids is explored by means of generalised (drift) fluid simulations. Recent 
numerical results on turbulence in dusty magnetised plasmas, strongly coupled fluids, semi-classical 
("quantum") plasmas and in rotating quantum fluids are reviewed and discussed. 
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I. INTRODUCTION 



Complex plasmas are characterised by strong coupling or complex interactions among 
the constituing particles (electrons, one or more ion species, and possibly neutral atoms and 
molecules), in particular in space and laboratory plasmas which are cooL dense, and/or may 
contain a significant amount of impurities, like massive dust particles 

Strongly coupled plasmas are specifically characterized by a larger than unity coupling 
coefficient, defined as the ratio between Coulomb energy and kinetic energy, and for example 
occur in laser-produced plasmas or in compact astrophysical objects, which in addition to 
strong coupling may for extreme parameters also manifest quantum effects. Dusty plasmas 
occur in interplanetary and interstellar space, in the edge of fusion experiments, and are 
studied by dedicated laboratory experiments I2I. An other example for complex plasmas 
of topical interest are ultra-cold neutral plasmas , obtained in the laboratory by laser 
ionization of ultra-cold matter like Bose-Einstein condensates. 

The physics of complex plasmas has received considerable attention over recent years. 
However, the theory and simulation of nonlinear dynamics and turbulent structures in mag- 
netised complex plasmas so far remains rather unexplored - despite the fact that static or 
dynamic magnetic fields are ubiquitously present in many space and laboratory plasmas and 
lead to an abundance of additional dynamical features. 

Then again, turbulence, vortices and fiows in magnetised ideal (high temperature) plas- 
mas have always been a topic of considerable interest for fusion plasma physics in general, 
and in particular for the fusion related previous work of the author of this contribution jsl. 
In this contribution we are going to explore the formation and evolution of nonlinear and 
turbulent dynamical structures in magnetised complex plasmas by means of generalised drift 
fiuid simulations. 

It should be stated clearly that hydrodynamic models for plasma turbulence, and specif- 
ically in complex plasmas, most often can only be considered as a first approximation to 
the problem. The importance of kinetic effects is rather a rule than an exception, and for 
magnetised plasmas a 5-d gyrokinetic model can in many cases be regarded as more appro- 
priate, in particular if the coUisionality is low, particle trapping becomes important, or other 
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kinetic effects enter into the picture. If the relevant mode frequencies are not much lower 
than the ion gyrofrequency then also the gyrokinetic model would have to be replaced by a 
full Vlasov-Maxwell or N-particle model. These are in general computationally very inten- 
sive and often not affordable for many realistic problem scales. With noteable exceptions: 
recent ultra-cold neutral plasma experients, for example, consist usually only of a number 
of particles in the order of 10^ which is just manageable by advanced particle codes. For 
quantum condensate matter the numerical solution of the Gross-Pitaeveskii equation (GPE) 
(based on the nonlinear Schroedinger equation) is standard and delivers the full quantum 
physics involved. Such more realistic models, if available, should not be given up without 
need on cost of simplified fluid models. 

Then again, simplified hydrodynamic models can be quite instructive in terms of didac- 
tical aspects and may lend a more intuitive approach to the problem. It is worthwile to 
investigate, in which minimal models certain effects still qualitatively appear (even if the 
results may be quantitatively differing from the fundamental models). And moreover, for 
some situations more fundamental models than the fluid picture may even not be available 
(yet), or, as noted above, too expensive to be solved numerically. The bottom line is that 
fluid models may be instructive, but the user should be aware about their limitations. 

Along this motivation, the Hasegawa-Wakatani (HW) model for drift wave turbulence 
in magnetised plasmas is still often refered to for didactical purposes, even since more so- 
phisticates (gyrokinetic or gyrofiuid) models have become generally available. HW is the 
minimal model allowing for fundamental insights into the drift instability mechanism, non- 
linear vortex development, and emergence of zonal flows out of drift wave turbulence. The 
numerical solution of the 2-d HW equations is cheap and feasible even on a contemporary 
laptop. We here want to use modified versions of the HW equations for a first a proach to 
turbulence in complex magnetised plasmas. 
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II. 2-DIMENSIONAL TURBULENCE MODEL 

The HW model for resistive drift wave turbulence js] accounts for nonlinear instability 
driven by a gradient Vno(x) in plasma density and resistive parallel coupling between fluc- 
tuations of density n and electrostatic potential 0. The resulting turbulent state of ExB 
vortices in the (x, y) drift plane perpendicular to the magnetic field can form low-frequency 
ky = zonal flow structures with a finite wave number k^. 

The model is usually derived from isothermal electrostatic two-fluid equations for elec- 
trons and ions. Under drift approximation the perpendicular momentum equations deliver 
the low-frequency fluid drift velocities, notably the ExB drift, diamagnetic drift and polar- 
isation drift. In first order the ExB drift velocity enters into the nonlinear advection term 
of the density (continuum) equations, and the polarisation drifts enters through its finite 
divergence. 

Here we want to sketch out another approach and briefly motivate the derivation of the 
HW equations from gyrokinetic theory. The gyrokinetic equation is an evolution equation 
for the 5-d distribution function f{-x,v\\,v±,t) with respect to guiding center coordinates in 
a magnetised plasma. Modern gyrokinetic theory 7| is based on first stating the problem in 
Hamiltonian formalism and Lie transforming the corresponding Euler-Lagrange equations 
to eliminate the gyroangle coordinate under the assumption that the gyromotion is fast 
compared to all other relevant time scales {u <^ Qi). The moment expansion into (gyro)fluid 
equations based on this procedure involving a gyrokinetic Hamiltonian ideally conserves 
energy. 

Drift ordering further enters via smallness of the drift scale ps <^ Lj_ compared to back- 
ground gradient lengths and an ordering ^ A;|| with respect to orientation perpendicular 
and along the magnetic field. A minimal form of a gyrokinetic equation may be constructed 
by neglecting finite Larmor radius (FLR) effects, neglecting parallel dynamics, assuming a 
homogeneous magnetic field, and keeping only ExB dynamics in the convection, to obtain 
dtf + y^EccB ■ V/ = 0. 

Integration over velocity space is trivial when the most simple (fluid) model for the dis- 
tribution function / = uFq as a product between macroscopic density n and a Maxwellian 
Fq is assumed. This results in advection equations dtUs + ^exB ■ = for the densi- 
ties Ug of species s = {e,i). These are coupled by the polarisation equation — riiC = 
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{norrii / B'^)'V\(f) = fl, relating the densities to the vorticity Q. 

The ion density equation is usually replaced by a vorticity equation, which is obtained 
by subtracting the density equations. For the parallel dynamics again the simplest possible 
model is assumed, relating the electron current for cold ions electrostatically to the density 
and potential via Ohm's law: J\\ = (l/r7||)[(Te/nee)V||ne — 

The resulting quasi-two-dimensional equations are normalised as 

5oe0/Te — > 0, 6one/no ^ n + N{x), tCs/L±^t, x/L±^x, (1) 

with 60 = Ps/L±, cj. = Te/rrii, and = rriiT^/ {eBY . The density has been split into a 
fluctuating component n and a constant background N{x). The dissipative coupling through 
the current is parametrised by -D = k'^^{L^/ Cg) / {T^/noe^rj\\). The resulting set of equations 
is the HW standard model for resistive drift wave turbulence: 



dtn + [0,n] = -[(t),N]+ D{(t)-n) (2) 

dtn + [(p,n] = D{(p-n) (3) 

with Vi0 = Vl. (4) 
The general properties of the HW model have been extensively discussed elsewhere (e.g. 



in Refs. |8l-ll8|). The HW model supports unstable drift waves and saturated turbulence by 
tapping the free energy in the background density gradient N{x) through resisive coupling 
via D. 

The hydrodynamic limit of the two-dimensional Euler equation is recovered for D = 
and = const, while the adiabatic limit D ^ 1 asymptotically corresponds to the 
Hasegawa-Mima-Charney-Obukhov equation. The Hasegawa-Mima (HM) system 19|] in 
itself is stable and may be used to study decaying turbulence, or it could be artificially 
driven. HM is isomorphic to the Charney-Obukhov equations for rotating 2-d fluids (like 



planetary atmospheres) including Rossby wave dynamics 20|. 

A characteristic property of these two-dimensional fluid systems is the possibility for 
brmation of large-scale zonal flow structures that are coupled to the turbulent spectrum 



2l| . which is a manifestation of the dual cascade nature of 2-d turbulence. 

The importance of this set of equations for plasma physics has been underlined by the 
European Physical Society in awarding the 2011 Hannes Alfven prize to Hasegawa, Mima 



and Diamond. Which is of course a bit ironic, as the HW/HM models exphcitly neglect 
Alfven dynamics. 

We numerically solve equations ([2]) and with an explicit 3rd order Karniadakis time 



stepping scheme 



22| . and the Poisson brackets [a, b] = (dra)(d yb) — {dya){djj) are evaluated 



with the energy and enstrophy conserving Arakawa method 23|. The numerical method 



is equivalent to the one introduced in Refs. 



24 



25|. Hyperviscuous operators v^V^, with 



z/"^ = —2 ■ 10~^, are added for numerical stability to the right hand side of both equations 
([2]) and ([3]), acting on n and f2, respectively. We solve eq. (jl]) in A; space by evaluation of 
(pk = — ^fc/^i employing the FFTW3 transform. The equations are discretised on a 2-d 
rectangular (x, y) grid with various (in general not quadratic) box dimensions. Boundary 
conditions are periodic in y and either periodic or Dirichlet in x. An example for a typical 
vorticity field VL{x,y) in fully developed drift wave turbulence is shown figure 1. 



. V - 



*■ .A 







FIG. 1: Yovticity 0,{x,y) in fully developed Hasegawa-Wakatani drift wave turbulence. Negative 
values are colored in blue, positive values in red. The physical domain is 64ps x 128/?s with a 
resolution of Ux = 512 and Uy = 1024. The figure is here rotated by 90 degrees for better display. 
The vortices can often be identified to appear in dipolar pairs. 
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III. 2-D TURBULENCE IN COMPLEX PLASMAS 



The motivation for the present contribution originated from discussions during a visit 
of Padma K. Shukla to University of Innsbruck in 2010. It became clear that a number of 
problems regarding waves and instabilities in complex plasmas, that Shukla and many other 
authors had recently addressed by means of linear theory or in 1-d, could often straightfor- 
wardly be generalised to a 2-d nonlinear HW like system of fluid equations and put into a 
form which is directly treatable by our numerical scheme. Shukla had in particular suggested 
to investigate turbulence in 2-d models for quantum magnetoplasmas, dusty plasmas and 
strongly coupled fluids. In the following we review the initial results of our simulations. 



A. Drift wave turbulence in the presence of a dust density gradient 

I 

Shukla and Varma [261] have proposed in 1993 a model to study the effect of static, 
immobile dust grains on waves and instabilities in plasmas. It had been shown that the 
presence of a static dust density gradient results in modes with a frequency proportional to 
the dust density gradient scale. 

This model has now been cast into a 2-d HW like form that allows the study of drift wave 
turbulence in the presence of a density gradient of immobile charged dust particles: 

dtn + [(j),n] = -{b + a) dy(f) + D{(p-n), (5) 
dtn + [0, n] = -ady<p--iVL + D{<p - n), (6) 

where a = {epdZdc/ B()nio)dxndo with e = (+1, —1) for (negative, positive) dust, pd = Cd/^d, 



where q = {akBTe/rriiY/'^ is the modified ion-sound speed 27| with a = njo/''^eo > 1- The 
dust related viscosity is given by 7. The constant plasma density gradient (derived from 
N{x)) is here parametrised by 6 = L±/Ln, which in our present normalisation is unity. 
This set of modified HW equations has been derived and numerically solved by Shukla and 



Kendl in ref. [28[. The quasi-linear drift wave instability through an ExB vortex growing 
out of an initial Gaussian sensity perturbation has been studied for typical resistive drift 
wave parameters and various dust gradient scales a. It has been found that the presence of a 
co-aligned density gradient of positively charged dust (a > 0) strongly enhances the resistive 
drift wave instability, while a counter-aligned dust gradient leads to a damping of the drift 
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waves. The same effect has been seen on fully developed turbulence: co-aligned gradients 
result in larger fluctuation amplitudes and more pronounced small-scale structures at high 
k, whereas counter-alignment damps the turbulence into quasi-linear modes 



The model so far has been restricted on immobile dust grains and did not consider 
dust charging effects. A comprehensive 2-d multi-species gyrofluid code for a more general 
modelling of dusty plasma turbulence is presently being developed by the author. 

B. Generalised (viscoelastic) hydrodynamics 

Strongly coupled fluids are often modelled in "generalised hydrodynamics" by including 
a viscoleastic relaxation time r as lowest order manifestation of kinetic coupling effects in 
the (incompressible) Navier-Stokes equation: 

(1 + r dt) dtVl = aV\n (7) 

For numerical treatment with our existing solvers we split this equation in two dimensions 
into the coupled set 

dtn+[<P,n] = -[l/r){^l-^) (8) 

dt^ + [0, n] = avln (9) 

so that for r — )• the Navier-Stokes equation is asymptotically recovered. We have stud- 
ied vortex stability, propagation and decay with this system. Remarkably, we have found 
for comparable parameters qualitatively nearly identical results as in a recent study ^] 
that employed a "first principles" molecular dynamics code to simulate coherent vortices in 
strongly coupled liquids. 

A personal favourite is the nonlinear evolution of a vortex with just the right relation 
between initial amplitude and radius (i.e. half- width of a Gaussian). It is well known from 
classical hydrodynamics that a strong vortex can develop into a rotating tripole and further 



perform along various routes until decay 



30|. 



For certain parameters, the tripole can split up into a dipole pair and a single monopolar 
vortex. For some parameters the single may just hastily exit the system and leave the pair, 
which is further waltzing around each other alone. 

For other initial parameters, the three can perform a complicated kind of polka, where 
the separated monopole circles around the rotating dipole in some distance, approaches the 
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pair again to exchange partners, and the freely relased vortex is now circhng around the 
new pair... and so on to viscosity. 

A novel aspect regarding this dance, which we now studied with the above set of equations, 
is the combined effect of viscosity a and viscoelastic relaxation time r on this triple polka. 
We find that by adding high viscosity (e.g. a = 0.1) the initial vortex may rather (again 
depending on initial amplitude and scale) disperse outwards and not further develop any 
nonlinear structures. Adding finite r (here in the range between 1 and 10) can compensate 
for the dispersive action of a, and the dipolar or tripolar vortex structures of the inviscid 
case can be recovered. 



C. De Broglie screening effect in semi-classical plasmas 



Recently, Shukla and Kendl have derived and numerically solved a semi-classical generali- 



sation of the E 
corrections 



W equations for dense degenerate Fermi plasmas including quantum pressure 



3l|. The model is based on the quantum magneto-hydrodynamic model |32l. |33|. 



The dispersion relation of drift waves in magnetised quantum hydrodynamic plasmas has 



before been studied by a number of authors in similar linear models, e.g. refs . | 34 
The semi-classical HW equations in an inhomogeneous magnetic field are 

dtVL + [0, fi] = L)(0 - An) - K{K*n) 
dtn + [0, n] = D{(j)- An) + K{(f) - A*n) 



(10) 
(11) 



where the curvature operator is defined as 

Kif) = -cV X (B/i?2) . V/ = -cV X [(B X Vf)/B'] . 

The quantum corrections effect enters through A = 1 — /3^V^ and A* = 1 -|- /3^V^, with 
/3 = \q/ps, where A, = h/ y/^m^Tp is the electron de Broglie length and Ps the drift scale, 
both here defined at Fermi temperature. 

The novel feature of this system is that we find a finite de Broglie length (FBL) screening 
effect on density fluctuations, which is analogous to the well-known FLR gyro screening and 
Debye screening of plasma turbulence. 

An important caveat concerns the range of validity of the quantum hydrodynamic model. 
The fluid-like equations including the quantum pressure are derived from the nonlinear 



Schroedinger equations by means of an eikonal ansatz, replacing the quantum mechanical 
wave function by an amplitude n ~ which follows a continuity equation, and a phase 
function S", which dynamically evolves in space and time like a fluid velocity in Euler's 
equation. It is important to remember that the eikonal ansatz is valid only in the long wave 
length regime where all wave lengths 271 /k ^ Ag are required to be much larger than the 
de Broglie length. This is analogous to the relation between electromagnetic wave optics 
and geometric ray optics. While this remark may perhaps seem obvious to most readers, 
it has to be noted that some authors applying quantum plasma hydrodynamic models do 
unfortunately not seem to be aware of these restrictions. 

For quasi-linear drift waves the relevant wave lengths are in the order of a few drift scales 
Psi with most unstable modes usually found around 0.1 < Psk < 1. Choosing the parameter 
/3 = \q/ ps well below unity is therefore a safe bet for linear theory. In fully developed 
turbulence the cascade however in principle reaches down to viscous scales. We do not have 
a good model of where to properly insert a viscous or Reynolds cut-off in dense quantum 
plasma turbulence. Cutting the spectrum off by hyperviscosity at scales related to ps/4 
seems tolerable in our experience. In any case we require a small quantum parameter /3 <C 1 
for consistency. 

In our recent quantum effect studies on resistive drift wave turbulence 3l|] and interchange 



modes |36l we have used /3 = 0.5 as a (perhaps already somewhat questionable) upper limit 
which is barely consistent with the model validity. The main message here is that our results 
are highly interesting from an academic point of view, regarding the newly discovered FBL 
screening effect. From a practical perspective, the quantum parameter in a valid range /3 ^ 1 
gives only minor corrections to the classical results. When compared to gyrofiuids, this limit 
would somewhat correspond to a Taylor approximation of the gyro- averaging operators - 
which also is either practically insignificant for most purposes if used in the correct limit, or 
otherwise simply gives a wrong result. 

The validity of quantum hydrodynamics starts to cease at scales just where it is beginning 
to become interesting. In particular, the important phenomenon of vortex quantisation is 
not self-consistently treatable at all with the hydrodynamic equations (other than being 
artifically modelled as specific point or disc vortices). 
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D. Quantum Charney-Obukhov turbulence 



The same notes of caution apply to the last quasi- 2-d turbulent system, which we here 
present and have recently addressed numerically. As mentioned in the introduction, there 
is a well-known isomorphism between the Hasegawa-Mima and Charney-Obukhov models. 
The latter describes rotating fluids in the presence of a gradient in rotation frequency, which 
is applicable as a basic model for atmospherical dynamics, and as solution includes Rossby 
waves and the secondary formation of jet streams. Recent experiments on strongly rotating 
Bose-Einstein condensates provided a motivation to apply the quantum hydrodynamic ap- 
proach including strong rotation. A quantum Charney-Obukhov equation has recently been 
suggested and linearly analysed in ref. js?]: 



where i/j = A(j) with A = 1 — /3^V^. This equation can be treated with our numerical 
scheme by explicitly evaluating dW = [A0, + V^A0 — InA^], and inverting Wk = (1 + 
k"^ + (l/2)/3^A;^)0fc in Fourier space. We have numerically solved this system by using a 
Gaussian background rotation N{r) to model a rotating condensate in a trap and initially 
superimposing a shear flow perturbation at half radius. The shear flow is found to be 
breaking up into several vortex-like modes rotating around the center of the condensate 
(see figure 2), and then is rather rapidly decaying, when the shear flow perturbation is not 
continously driven. 

To the best of our knowledge, such a nonlinear 2-d quantum hydrodynamical model of 
rapidly rotating Bose-Einstein condensates has not been simulated before. But again, in 
the valid range the small quantum factor /3 does not produce any significant effects. The 
behaviour of the system is found to be very similar to a classical rotating fluid. 

What of course is completely absent in this model is vortex quantisation. And as this 
is just the whole crux of quantum turbulence, the hydrodynamic approach on BECs can in 
the personal opinion of the author be regarded as rather useless. Nothing new is learned, 
much is lost. There seems to be no good reason for not rather solving the CPE (which is 
really rather standard now) to model BECs. 




(12) 
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FIG. 2: A shear How perturbation is superimposed on a circular rotating quantum Quid in a 2-d 
numerical solution of the (actually semi-classical) quantum Charney-Obukhov model. Here the 
perturbed stream function (p{x, y) is shown, which develops a modal structure and (depending on 
parameters, here (3 = 0.8) either rapidly decays or transiently forms a chain of rotating vortices. 
Note the absence of vortex quantisation in this simplified model. The color coding is blue for 
positive (corresponding to the initial perturbation) and red for negative values of (p. Stream lines 
of the Euid How are on isocontours of the stream function. 



IV. CONCLUSIONS 

We have numerically studied the formation and evolution of nonlinear and turbulent 
dynamical structures in two-dimensional complex plasmas and fluids. Generalised (drift) 
fluid simulations have been developed on the basis of the Hasegawa-Wakatani model and 
the 2-d Navier-Stokes vorticity equation. 

Recently published results on turbulence in dusty magnetised plasmas and semi-classical 
("quantum") plasmas have been reviewed. 

New results on 2-d generalised hydrodynamics including viscoelastic relaxation effects 
were presented. The first nonlinear quantum hydrodynamic simulation of a rotating Bose- 
Einstein condensate has been presented. The limits of validity of the quantum hydrodynamic 
model with respect to turbulence was critically discussed. 
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